clear all

% Simulate (gamma = 2 x estimated value), partial equilibrium

global Dividend

Dividend=0;

% Given Par2 optimal:
n = 50; tau = 0.1;
r = 0.036;
Par1 = zeros(25,1);
load ParEst_eh0_time0_Mex ParEst_eh0_time0_Mex
load ParEst_eh0_time1_Mex ParEst_eh0_time1_Mex
Par1(1:24) = ParEst_eh0_time0_Mex;
Par1(21) = 10^Par1(21);
Par1(22:24) = Par1(22:24)*1000;
Par1(25) = 2*500*ParEst_eh0_time1_Mex(1); 
gamma =  Par1(25); a = Par1(22); b1 = Par1(23); b2 = Par1(24);
alphaf = Par1(17); betaf = Par1(18); alphai = Par1(19); betai = Par1(20); theta = Par1(21);
load W_time1_age0_eh0_comp11.out
w_f_g=W_time1_age0_eh0_comp11; % Read wages, g and f distributions
wf_1=w_f_g(:,1);        % Wage in formal sector
wi_1=w_f_g(:,2);        % Wage in informal sector
wf_2=w_f_g(:,7);      % Wage in formal sector
wi_2=w_f_g(:,8);      % Wage in informal sector

[Ff_1,ff_1,Fi_1,fi_1,Ff_2,ff_2,Fi_2,fi_2]=solve_wage_offer_dist(Par1,wf_1,wi_1,wf_2,wi_2,n);
[Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn,Rw_fi_ni, Rw_ii_ni, Rw_ii_in, Rw_if_in,...
     Rw_fi_fn, Rw_ff_fn, Rw_if_nf, Rw_ff_nf, Rw_in_nn, Rw_fn_nn, Rw_ni_nn, Rw_nf_nn]=solve_value_function_A3_C(Par1,Ff_1,wf_1,Fi_1,wi_1,Ff_2,wf_2,Fi_2,wi_2,n);
[HHstates,wage_moments,transitions,gff, gfi, gif, gii, gfn, gnf, gin, gni]=data_simulation_A12_end(Par1,ff_1,wf_1,fi_1,wi_1,ff_2,wf_2,fi_2,wi_2,Wff,Wfi,Wfn,Wif,Wnf,Wii,Win,Wni,Wnn);
wage_moments_gamma_1_PE = wage_moments;
Rw_1_f_gamma_1_PE = Rw_fi_ni; Rw_1_i_gamma_1_PE = Rw_ii_ni;  Rw_2_i_gamma_1_PE = Rw_ii_in; Rw_2_f_gamma_1_PE = Rw_if_in;

save Rw_1_f_gamma_1_PE Rw_1_f_gamma_1_PE
save Rw_2_f_gamma_1_PE Rw_2_f_gamma_1_PE
save Rw_1_i_gamma_1_PE Rw_1_i_gamma_1_PE
save Rw_2_i_gamma_1_PE Rw_2_i_gamma_1_PE
HH_states_gamma_1_PE = HHstates;
save wage_moments_gamma_1_PE wage_moments_gamma_1_PE
save HH_states_gamma_1_PE HH_states_gamma_1_PE

mff=HHstates(1);
mfi=HHstates(2);
mfn=HHstates(3);
mif=HHstates(4);
mnf=HHstates(5);
mii=HHstates(6);
min_=HHstates(7);
mni=HHstates(8);
mnn=HHstates(9);
Un_Head = mnf + mni + mnn; Un_Spouse = mnn + min_ + mfn; Inf_rate = (mnn + min_ + mni + mii);
Agg_stocks_gamma_1_PE= [Inf_rate; Un_Head; Un_Spouse];
save Agg_stocks_gamma_1_PE Agg_stocks_gamma_1_PE
gf_1 = zeros (50,1); gf_2 = zeros (50,1); gi_1 =  zeros (50,1); gi_2 = zeros (50,1);
 for i = 1:50
 gf_1(i) =  (mff*(sum(gff(i,:))) + mfi*(sum(gfi(i,:))) + mfn*gfn(i))/(mff+mfi+mfn);   
 end

  for i = 1:50
 gi_1(i) =  (mif*(sum(gif(i,:))) + mii*(sum(gii(i,:))) + min_*gin(i))/(mif+mii+min_);   
  end
 
   for i = 1:50
 gf_2(i) =  (mff*(sum(gff(:,i))) + mif*(sum(gif(:,i))) + mnf*gnf(i))/(mff+mif+mnf);   
   end
 
    for i = 1:50
 gi_2(i) =  (mfi*(sum(gfi(:,i))) + mii*(sum(gii(:,i))) + mni*gni(i))/(mfi+mii + mni);   
 end

% have to construct the joint probabilities (for ex, gff(wf_1,wf_2)) or integrate over the
% marginal distributions gf_1 and gf_2, whichever is easier..

Welfare_total=r*(sum(sum(mff*Wff.*gff+mfi*Wfi.*gfi+mif*Wif.*gif+mii*Wii.*gii)))+r*sum(mfn*Wfn.*gfn+mnf*Wnf.*gnf+min_*Win.*gin+mni*Wni.*gni)+r*mnn*Wnn;
Welfare_formal1=(r*(sum(sum(mff*Wff.*gff+mfi*Wfi.*gfi)))+r*sum(mfn*Wfn.*gfn))/(mff+mfi+mfn);
Welfare_informal1=(r*(sum(sum(mif*Wif.*gif+mii*Wii.*gii)))+r*sum(min_*Win.*gin))/(mif+mii+min_);
Welfare_nonemployed1=(r*sum(mnf*Wnf.*gnf+mni*Wni.*gni)+r*mnn*Wnn)/(mnf+mni+mnn);
Welfare_formal2=(r*(sum(sum(mff*Wff.*gff+mif*Wif.*gif)))+r*sum(mnf*Wnf.*gnf))/(mnf+mif+mff);
Welfare_informal2=(r*(sum(sum(mfi*Wfi.*gfi+mii*Wii.*gii)))+r*sum(mni*Wni.*gni))/(mni+mii+mfi);
Welfare_nonemployed2=(r*sum(mfn*Wfn.*gfn+min_*Win.*gin)+r*mnn*Wnn)/(min_+mfn+mnn);
Welfare_gamma_1_PE=[Welfare_total;Welfare_formal1;Welfare_informal1;Welfare_nonemployed1;Welfare_formal2;Welfare_informal2;Welfare_nonemployed2]; % IN THOUSANDS, IF NORMALIZED
save Welfare_gamma_1_PE Welfare_gamma_1_PE

% Generate LOG WAGES by the F distribution 
[wvector_f_1_est,Fvector_f_1_w]=main_quantiles(1-Ff_1,wf_1,n);
[wvector_i_1_est,Fvector_i_1_w]=main_quantiles(1-Fi_1,wi_1,n);
[wvector_f_2_est,Fvector_f_2_w]=main_quantiles(1-Ff_2,wf_2,n);
[wvector_i_2_est,Fvector_i_2_w]=main_quantiles(1-Fi_2,wi_2,n);
mean_wf_1f=log(sum(wf_1.*ff_1));
mean_wi_1f=log(sum(wi_1.*fi_1));
mean_wf_2f=log(sum(wf_2.*ff_2));
mean_wi_2f=log(sum(wi_2.*fi_2));
wage_moments_under_F=[[wvector_f_1_est;mean_wf_1f],[wvector_i_1_est;mean_wi_1f],[wvector_f_2_est;mean_wf_2f],[wvector_i_2_est;mean_wi_2f]];
save wage_moments_under_f_gamma_1_PE wage_moments_under_F

% Transitions 
Trans_head = [transitions(3); transitions(4); transitions(1); transitions(5); transitions(2); transitions(6); transitions(13); transitions(14)];
Trans_spouse = [transitions(9); transitions(10); transitions(7); transitions(11); transitions(8); transitions(12); transitions(15); transitions(16)];
Transitions_gamma_1_PE = [Trans_head,Trans_spouse];
save Transitions_gamma_1_PE Transitions_gamma_1_PE

% Constructing produtivities for simulations

% marginal distributions
mf_1=mff+mfi+mfn;
mi_1=mif+mii+min_;
mf_2=mff+mif+mnf;
mi_2=mfi+mii+mni;


% firm size
el_f_1=mf_1*gf_1./ff_1; 
% el_f_1(1)=el_f_1(2); % may need smoothing
% el_f_1(n)=el_f_1(n-1);
el_i_1=mi_1*gi_1./fi_1;
% el_i_1(1)=el_i_1(2); % may need smoothing
% el_i_1(n)=el_i_1(n-1);
el_f_2=mf_2*gf_2./ff_2; 
% el_f_2(1)=el_f_2(2); % may need smoothing
% el_f_2(n)=el_f_2(n-1);
el_i_2=mi_2*gi_2./fi_2;
% el_i_2(1)=el_i_2(2); % may need smoothing
% el_i_2(n)=el_i_2(n-1);

load fitted_model p_f_1 p_i_1 p_f_2 p_i_2 Ff_1 Fi_1 Ff_2 Fi_2 ff_1 fi_1 ff_2 fi_2

% M (total employed workforce) in 2005 = 41441076
% Fraction of firms in the formal sector in 2002 includes self-emp. (ENAMIN/INEGI: 23.7%)
% Number of formal firms (IMSS registry, BOSCH AND CAMPOS-VASQUEZ 2005): 800,000
% Number of firms = 800000/0.237
M_N=41441076/(800000/0.237);

pi_f_1=(p_f_1-wf_1*(1+tau)).*el_f_1*M_N/0.237;
pi_i_1=(p_i_1-wi_1).*el_i_1*M_N/0.763;
pi_f_2=(p_f_2-wf_2*(1+tau)).*el_f_2*M_N/0.237;
pi_i_2=(p_i_2-wi_2).*el_i_2*M_N/0.763;
% may need smoothing
% pi_f_1s=smoothdata(pi_f_1(1:25),'rlowess',20);
% pi_i_1s=smoothdata(pi_i_1(1:18),'rlowess',20);
% pi_f_2s=smoothdata(pi_f_2(1:25),'rlowess',20);
% pi_i_2s=smoothdata(pi_i_2(1:18),'rlowess',20);

profitFormalFirms_1=1./(M_N*mf_1/0.237)*sum(pi_f_1(1:25).*ff_1(1:25)); % per worker 
profitInformalFirms_1=1./(M_N*mi_1/0.763)*sum(pi_i_1(1:18).*fi_1(1:18)); % per worker 
totalprofitFirms_1=mf_1*profitFormalFirms_1+mi_1*profitInformalFirms_1;  
profitFormalFirms_2=1./(M_N*mf_2/0.237)*sum(pi_f_2(1:25).*ff_2(1:25)); % per worker 
profitInformalFirms_2=1./(M_N*mi_2/0.763)*sum(pi_i_2(1:18).*fi_2(1:18)); % per worker 
totalprofitFirms_2=mf_2*profitFormalFirms_2+mi_2*profitInformalFirms_2;  
Profit_gamma_1_PE=[totalprofitFirms_1;profitFormalFirms_1;profitInformalFirms_1;NaN;totalprofitFirms_2;profitFormalFirms_2;profitInformalFirms_2]; % per worker
save Profit_gamma_1_PE Profit_gamma_1_PE

% normalized firm size
[elvector_f_1,Fvector_f_1]=main_quantiles2(1-Ff_1,el_f_1,n)
[elvector_i_1,Fvector_i_1]=main_quantiles2(1-Fi_1,el_i_1,n)
[elvector_f_2,Fvector_f_2]=main_quantiles2(1-Ff_2,el_f_2,n)
[elvector_i_2,Fvector_i_2]=main_quantiles2(1-Fi_2,el_i_2,n)
mean_el_f_1=sum(el_f_1.*ff_1);
mean_el_i_1=sum(el_i_1.*fi_1);
mean_el_f_2=sum(el_f_2.*ff_2);
mean_el_i_2=sum(el_i_2.*fi_2);
Firm_size_gamma_1_PE=[[elvector_f_1,elvector_i_1,elvector_f_2,elvector_i_2];[mean_el_f_1,mean_el_i_1,mean_el_f_2,mean_el_i_2]];
save Firm_size_gamma_1_PE Firm_size_gamma_1_PE

